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ABSTRACT 



Context. We present an improved version of the SimpleX method for radiative transfer on an unstructured Delaunay grid. The grid 
samples the medium through which photons are transported in an optimal way for fast radiative transfer calculations. 
Aims. We study the detailed working of SimpleX in test problems and show improvements over the previous version of the method. 
Methods. We have applied a direction conserving transport scheme that correctly transports photons in the optically thin regime, a 
regime where the original SimpleX algorithm lost its accuracy. In addition, a scheme of dynamic grid updates is employed to ensure 
correct photon transport when the optical depth changes during a simulation. For the application to large data sets, the method is 
parallellised for distributed memory machines using MPI. 

Results. To test the new method, we have performed standard tests for cosmological radiative transfer. We show excellent correspon- 
dence with both the analytical solution (when available) and to the results of other codes compared to the former version of SimpleX, 
without sacrificing the benefits of the high computational speed of the method. 
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1. Introduction 

In many astrophysical applications radiative processes play an 
important, and occasionally dominant, role. The interaction be- 
tween radiation and matter leads to possibly strong feedback ef- 
fects on the medium caused by heating and cooling, radiation 
pressure and the change of the ionisation and excitation state of 
the medium. It is therefore of crucial importance to include these 
effects in the simulations of hydrodynamic flows. However, ra- 
diative transfer is a very complex process due to the high di- 
mensionality of the problem. The specific intensity /(O, x, t, v) 
depends on seven variables, and is impossible to solve for in an 
analytic way in general. With the exception of certain limiting 
cases where an analytical solution exists, one therefore has to 
rely on numerical methods to obtain the desired solution. 

Several kinds of methods exist for this purpose, all of which 
have their specific advantages and shortcomings. A relatively 
straightforward way of solving the r adiative transfer equ a tion is 
. the long characteristics method (e.g. lMihalas & Mihalasl (I1984) 
), where rays connect each source cell with all other cells and 
\ the transfer equation is solved along every ray. Although this 
method is relatively precise, it is computationally demanding, 
requiring N 2 interactions for N cells. A solution to this un- 
fortunate scalin g is implemented i n the short characteristics 
method (e.g. iKunasz & Aued (Il988l) ), where only neighbour- 
ing cells are connected by rays. Not only does this make the 
method computationally less expensive than the long character- 
istics method, it is also easier to parallellise. A drawback of this 
method is that it is known to cause numerical diffusion in the 
solution. In recent years, hybrid schemes combi ning the ben- 
efits of both approache s have been developed (iRijkhorst et alJ 
l2006l:lTrac & Cenll2007b . Instead of direct integration along the 
characteristics, Monte Carlo methods use a stochastic approach 
by sending photon packets along the rays. The properties of the 
photon packets and their interaction with the medium are de- 



termined by sampling the distribution functions of the relevant 
processes. Moment methods solve the moments of the radiative 
transfer equation, which allows for a computational speed-up in 
certain opacity regimes. In these methods there is a trade-off be- 
tween computation time and numerical diffusion in the solution, 
depending on what method is used to obtain the closure relation. 

What almost all of these methods have in common is that 
they use a predefined grid on which the radiative transfer calcu- 
lation is done. Most often this grid is given by a hydrodynamics 
simulation, which is either a regular, adaptive mesh refinement 
(AMR) or smoothed particle hydrodynamics (SPH) grid. These 
grids are not optimised for radiative transfer but for hydrody- 
namics calculations, possibly resulting in unphysical behaviour 
of the numerical solution. Moreover, the computation time of 
almost all methods scales linearly with the number of sources, 
which severely limits the range of applications. Exceptions are 
moment methods that generally do not scale with the number of 
sources, but sacrifice prec ision by introducing numerical diffu- 
sio n (e.g. iGnedin & Abell d2001h : ICenl (l2002h ) and the method 
by lPawlik & Schayel ([2008), where a source merging procedure 
is used to avoid linear scaling with the number of sources. For 
many applications, the linear scaling of the computation time 
with the number of sources becomes prohibitive, for example 
when simulating scattering processes, where effectively every 
cell might become a source. In the case of simulations of the 
epoch of reionisation, which is the topic of the second half of 
this paper, it is necessary to include many sources. It is there- 
fore essential to use a method for which the computation time is 
independent of the number of sources. 

The approach to solve the radiative transfer equation taken 
in this paper is radically different from the methods described 
earlier. Instead of using a predefined grid, the SimpleX method 
calculates the grid for the radiative transfer problem from the 
properties of the physical medium through which photons will 
be transported. This leads to a computationally fast method that 
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does not scale with the number of sources, making it an ideal tool 
for simulations of the epoch of reioni sation. A previous versio n 
of the method has been described in iRitzerveld & Ickel ( 2006), 
where the general idea of transport on random lattices is laid 
down, with a small section on the application to cosmological 
reionisation. The first comprehensive set of tests of the method 
were performed f or the Radiative Transfer Comparison Project 
(llliev et al .1120061) . In this project, the focus lay on comparing 
the performance of all participating codes in the test problems, 
making an in-depth analysis of the SimpleX results impossible. 

The aim of this paper is to describe the improvements of the 
SimpleX method since these previous two papers. Recently, we 
have pe rformed a detailed study of the error properties of the 
method dKruip et al.ll2Q09h (KPCI09), which has led to some es- 
sential improvements to the method. The two main problems in 
ballistic transport that were addressed in KPCI09, decollimation 
and deflection, are minimised both by using an alternative trans- 
port scheme in the opacity regime where these problems occur 
and by adapting the grid to changes in the opacity during a sim- 
ulation. In addition, the algorithm has been parallellised for dis- 
tributed memory. In this paper, we describe the working of the 
improved SimpleX method and provide a detailed analysis of the 
algorithm in test problems focusing on cosmological radiative 
transfer. 

The format of this paper is as follows. In Sect. [2) we give 
an overview of the SimpleX method and a description of the 
new features of the method. We then describe the parallellisation 
strategy in Sect.[3]and present the computational scaling proper- 
ties of the algorithm. In Sect. |4] we focus on the specific appli- 
cation of the SimpleX algorithm to the problem of cosmological 
radiative transfer, and describe test problems for this application 
in Sect. [3 Finally, we present a summary in Sect. [6j 



2. The SimpleX method 

In this section, we describe the basics of the SimpleX method and 
specifically the new features that were added to improve the per- 
formance in the lower opacity regimes. For the sake of clarity we 
repeat some essential i nformati on that was p resented earlier in 
Ritzerveld & Icke ( 2006) and Ritzerveld ( 2007), which is neces- 
sary to appreciate the new features of the method. We start with 
a description of how the unstructured grid is created and how 
to optimise it for radiative transfer calculations with SimpleX. 
We then proceed by describing the different ways of transport- 
ing photons on this grid, governed by the physical properties of 
the problem at hand. 



2.1. Grid calculation 

At the basis of the SimpleX method lies the unstructured grid on 
which the photons are transported. The grid adapts to the phys- 
ical properties of the medium through which the photons travel 
in such a way that more grid points are placed in regions with a 
higher opacity. The result is a higher resolution in places where 
it's needed most, there where the optical depth is highest. 



2.1.1. Point process 

The placement of the grid points is a stochastic process based on 
the Poisson process, which can be defined as follows. Suppose 
N(A) is the number of points in a non-empty subset A of the 
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Fig. 1. Two-dimensional example of a point distribution based 
on a homogeneous Poisson process (left) and based on a non- 
homogeneous Poisson process (right). 

volume S c R d , with d the dimension. Then the probability to 
that A contains x points is 



0> = P(N(A) = x) 



n p \A\e- n pM* 



0,1,2, 



(1) 



The only parameter in this process is the point intensity n p , 
which is a global constant. Every region in the volume has the 
same probability that points are placed there, which in our case 
corresponds to a constant opacity. An example of a homoge- 
neous Poisson process is shown on the left hand side of Fig. [T] 

To account for different opacity regimes inside the compu- 
tational volume, we use the non-homogeneous Poisson process, 
defined as 



P(N(A) = x) = 



n p (A)\A\e~ n ^ A)lAlx 



x = 0,1,2,... 



where 

n p (A) = I n p (x)dx. 
Ja 



(2) 



(3) 



The point intensity function n p (x) follows the opacity of the 
medium on global scale, while on local scale the point distri- 
bution retains the properties of the homogeneous Poisson dis- 
tribution. An example of a point distribution based on the non- 
homogeneous Poisson process is shown on the right hand side of 
Fig.[T] An alternative, possibly more physically intuitive, recipe 
for constructing the non-homogeneous Poisson process can be 
written as 



n p (x) = O * f(n(x)\ 



(4) 



that is, by defining the grid point distribution as a convolution 
of a homogeneous Poisson process and a function of the possi- 
bly inhomogeneous medium density distribution n(x). It is this 
recipe for the non-homogeneous Poisson process that we use to 
construct the SimpleX grid. Grid points are placed by randomly 
sampling the correlation function f(n(x)). We will discuss the 
exact shape of the correlation function in Sect. 12.1.31 

2.1.2. The Delaunay triangulation 

The grid points thus created form the nuclei of a Voron oi tessel- 
lation Given a set o f nuclei xu the Voronoi tessellation (iDirichletl 
1 1850H Voronoi 1908h is defined as V = {Q}, in which 



Q = \ye 



\xi-y\\ < \\xj-y\\ V Xj ± x t 



(5) 
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In other words, this means that every point inside a Voronoi cell 
is closer to the nucleus of that cell than to any other nucleus. 
By joining all nuclei that have a common facet (an ed ge in 2D, 
a wa ll in 3D), we create the Delaunay triangulation (iDelaunayl 
Thus, every nucleus is connected to its closest neigh- 
bours. A 2D example of a Voronoi tessellation and the corre- 
sponding Delaunay triangulation is shown in Fig.O 

The Delaunay triangulation consists of simplices that fill the 
entire domain. A simplex is the generalisation of a triangle in R d , 
so a triangle in R 2 and a tetrahedron in R 3 . In a valid Delaunay 
triangulation, all simplices obey the empty circumsphere crite- 
rion. The circumsphere of a simplex is the unique sphere that 
passes through each of the vertices that make up the simplex. In 
a valid Delaunay triangulation, no vertex exists inside this cir- 
cumsphere. 

For Voronoi tesselations and Delaunay triangulations that 
are constructed from a point process based on a homogenous 
Poisson process, so-called Poisson Delaunay triangulations, it is 
possible to derive some general properties relevant for our ra- 
diative ^ai^fe^jnethod. Tlie^resjilts were mainly derived by 
iMilesI (119701 119741) and lMolled (Il989h . Two important proper- 
ties for our purposes are the average number of neighbours of a 
vertex and the average distance between two connected vertices. 
The expectation value for the number of neighbours of a typical 
vertex in R 2 and R 3 is 

E 2D (£) = 6 (6) 
and 

48tt 2 

E 3D (£)= — +2*15.54. (7) 

The expectation value for the distance between two connected 
vertices in R 2 and R 3 is 

and 

E3D(L) = ^(i) 1/3 ^ 1/3 ^ L237 < /3 - (9) 

Note that these values are only exact for Delaunay triangula- 
tions constructed from a homogeneous Poisson process, while in 
SimpleX we use the non-homogeneous Poisson process to place 
the grid points. Except for regions in the domain with strong gra- 
dients in the point density, on local scale the point distribution re- 
sembles a homogeneous point distribution quite well. Therefore 
the properties of the Poisson Delaunay triangulation give a good 



qualitative idea of the properties of the grid on which we perform 
our radiative transfer calculations. 

SimpleX is set up in such a way that once the point distribu- 
tion is created according to Eq. ©, the Delaunay triangulation is 
calculated by an external software package. It is therefore possi- 
ble to use any package that suits the application at hand. For all 
simulations presented in this paper, the Delaunay triangulation is 
calculated using the QHull packag^J This is a software package 
written in C that is able to calculate the Delaunay triangulation, 
the surfaces and the volumes of the simplices in u p to 8 dimen- 
sion s. QHull is based on the Quickhull algorithm dBarber et al.l 
119951) . using the convex hull property of the Delaunay triangula- 
tion. QHull has the advantages that it computes the Delaunay tri- 
angulation in optimal time 0(NlogN), it is very stable against 
floating point round off errors in case points lie very close to 
each other and it is easy to implement as modular plugin routine. 
One of the drawbacks of QHull is that it triangulates the entire 
point set in one call, so it's impossible to add or delete points 
after the triangulation has been computed. This results in extra 
computational overhead in the grid dynamics scheme presented 
in Sect. 12.1.41 However, the computation time of the triangula- 
tion is small compared to the computation time of the radiative 
transfer (see also Fig. [6]), so in the present case the extra compu- 
tational overhead is acceptable. 

2.1 .3. The correlation function 

In the previous discussion we have not specified the exact shape 
of the correlation function f(n(x)) with which we sample the 
density distribution of the medium. In order for the grid to adapt 
to the properties of the medium, the correlation function should 
be a monotonically increasing function in n(x). Thus, the dis- 
tance between two connected vertices will be smaller in regions 
with high density. From basic transfer theory, we know that the 
local mean free path in a medium relates to the local medium 
density in the following way: 

m = (io) 

n(x)cr 

where cr is the total cross section, cr = o> consisting of dif- 
ferent interaction cross sections cry. If we compare this to the 
expectation value of the Delaunay edge length in Eq. ® and 
Eq. © it follows that if we choose the correlation function 
f(n(x)) to sample the d-th power of the local medium density, 
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e.g. f{x) oc x d , the local mean free path scales linearly with the 
expectation value of the Delaunay edge length via a constant c: 



(L D )ix) = cA(x). 



(11) 



Equation (ITTI) is a global relation with a global constant c, given 
by 



c = £(d,D,N)cr. 



(12) 



Here, D is the size of the computational domain and N the num- 
ber of vertices. 

This sampling recipe works very well in physical media 
where the density fluctuations are small. However, if the den- 
sity fluctuations are significant, sampling the density to the d-th 
power will favour the high density regions, resulting in a possi- 
bly severe undersampling of the low density regions. This under- 
sampling can have serious consequences for the radiative trans- 
fer calculation, for instance by causing preferential dire ctions 
that lead radiation around low density regions. See Sect. 15.3.21 
for an example of this effect. Moreover, KPCI09 showed that 
large gradients in the grid point distribution lead to systematic 
errors when transporting photons on this grid. 

We therefore need a sampling recipe that retains the advan- 
tages of the adaptive grid by keeping the dynamic range as large 
as possible, while maximising the minimum resolution of the 
grid. This is achieved by defining a reference density noix) above 
which the density is no longer sampled to the d-th power but a 
different power a. Both noix) and a depend on the properties of 
the medium that needs to be sampled. The two sampling recipes 
are smoothly joined by taking the harmonic mean, resulting in 
the following sampling function: 



/(*(*)) = 



nix) 
noix) 



-d 



nix) 
noix) 



(13) 



This sampling recipe favours low density regions by sampling 
those with a higher (d-th) power. 

Following KPCI09 we can define a sampling parameter Q n 
that is a measure of the gradients in the point distribution: 



Qn 



n p jx) 
Vn p ix)' 



(14) 



Similarly, we can define a measure of the gradients in the density 
distribution: 



Qy = 



nix) 
Vnix) ' 



(15) 



Using Eq. (TT3l) we can find the relationship between Q n and Q y : 





Q n = 1 .44 
»:"= 0.99 pt 



Q n = 1 44 
(*"= 0.75 a m 




Q n = 1 44 
a = 0.1 a |lm 



a = 0.99 a lim 



a = 0.1 aj 



= 12.0 
= 0.99 a |im 



= 12.0 
= 0.75 a |ln 



= 12.0 
= 0.1 a lim 



Fig. 3. Example of the influence of the sampling parameters Q n , 
a and noix) on the sampling of a 2D density field. In all fig- 
ures the number of grid points is 20000, sampling parameters 
are plotted such that Q n increases from top to bottom and a de- 
creases from left to right. The values for a are taken relative to 
the limiting value for which noix) = 0. From left to right a val- 
ues are 0.99a/; m , 0.75au m and 0Aau m . First row: Density field 
with a single density peak and a 1/r 2 profile with Q y = 1.44 and 
the sampling according to the square of the density. Clearly this 
sampling recipe leads to a bias towards the high density peak. 
Second row: sampling with Q n = 1.44 and au m = 1.0. Third row: 
sampling with Q n = 5.0 and au m - 0.29. Fourth row: sampling 
with Q n - 12.0 and au m = 0.12. Increasing Q n leads to smaller 
gradients in the point distribution, eventually leading to a ho- 
mogeneous point distribution. Decreasing a implies increasing 
noix). This leads to a flatter point distribution close to the den- 
sity peak but also to undersampling in the low density regions. 



Qn = Qy 



y «+y~ 



yiay~ a ~ l + dy- d ~ l Y 



(16) 



where we have defined y = nix) /noix). A choice for Q n sets 
the maximum gradient in the point distribution and thus the dy- 
namic range in the simulation. The higher Q n , the smaller the 
gradients in the point distribution. An example of a density field 
and point distribution with different values for Q n is shown in 
Fig. [3] This figure shows that increasing Q n indeed results in 
a smoother point distribution, with more points placed in the 
low density regions. A high Q n therefore makes it less likely 
that errors due to undersampling in low density regions occur. 
However, increasing Q n also decreases the dynamic range in the 
simulation, resulting in a possible undersampling of high den- 
sity peaks for high Q n values. For optimal sampling one should 



therefore choose the lowest Q n value for which numerical errors 
due to undersampling in low density regions are within a pre- 
defined tolerance set by the requirements of the simulation. For 
a more elaborate discussion of the Q n parameter, correspond- 
ing numerical artefacts and their solutions we refer the reader to 
KPCI09. 

Since Q y is a fixed property of the medium density, a certain 
value of Q n only exists for specific combinations of noix) and 
a. For every value of Q n there is a maximum a at which noix) 
goes to zero. The influence of different a values on the point 
distribution for fixed Q n is shown in Fig. [3] A lower value of 
a implies a higher value for noix). The result is that the high 
density peak is less pronounced, due to the lower value of a 
while at the same time less grid points are present in the lowest 
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density regions, due to the higher value of no(x). The reason for 
this is that in this example there is a gradient in the medium 
density everywhere, so a higher no(x) results in an emphasis on 
the regions where the density is close to this reference value. 
Therefore, the lowest density regions receive less points. It is 
therefore crucial to choose an a value that ensures that no strong 
density gradients exists at n(x) < no(x). 

By choosing the sampling parameter Q n and the reference 
density hq(x) we can create a grid where the dynamic range is 
maximal without causing numerical artefacts due to undersam- 
pling of the low density regions or large gradients in the point 
distribution. One has to be careful that by sampling different 
opacity regimes in a different way, the interaction coefficient of 
Eq. (TT2t is no longer a global constant, but we can take care of 
this by the way the interaction at each grid point is accounted 
for. This will be discussed in Sect. 12.21 

2.1 .4. A dynamic grid 

In the previous section we described how the SimpleX grid is 
created according to the properties of the medium. In this dis- 
cussion, it was assumed that the medium is static and does not 
change during the radiative transfer calculation. However, in re- 
ality the properties of the medium change continuously under 
influence of, for example, gravity and radiation. For this reason, 
the SimpleX grid should be updated every time step in case of 
full radiation hydrodynamics simulations. In this paper we do 
not consider the application of SimpleX to radiation hydrody- 
namics but instead focus on post-processing static density fields. 
We will discuss the application of SimpleX to radiation hydro- 
dynamics simulations in future work. 

Even though the gas density is assumed to be static, the prop- 
erties of the medium might still change during a radiative trans- 
fer calculation. For example, photo-ionisation lowers the optical 
depth for ionising radiation. Changes in the optical depth lead 
to a deviation from the recipe for grid point distribution of Eq. 
(TT31) . In other words, the grid is no longer an optimal representa- 
tion of the physical properties relevant for the radiative transfer 
calculation. KPCI09 showed that this might have serious conse- 
quences for the transport of photons through regions where the 
optical depth between grid points is much lower than unity. In 
Sect. 12.2.31 we describe a new transport scheme that minimises 
errors in the optically thin regime. This section describes a dif- 
ferent solution for transport through regions where the optical 
depth has severely changed, the removal of grid points. 

One of the advantages of the adaptive grid that SimpleX uses 
is that resolution is put where it is needed most, in the regions 
with highest optical depth. If during a radiative transfer calcula- 
tion the optical depth changes, we are effectively wasting com- 
putational resources in the regions with high resolution where 
the opacity has decreased. The high resolution is no longer nec- 
essary, since no interesting radiative transfer effects that need to 
be resolved at high resolution are taking place. The superfluous 
grid points only add to the computation time. We therefore re- 
move unnecessary grid points from the regions where the optical 
depth has significantly decreased. 

Another reason for the removal of superfluous grid points is 
that the transport of photons between grid points with an opti- 
cal depth lower than unity is prone to numerical errors. KPCI09 
showed that photons that travel through these regions are subject 
to decollimation and deflection. These effects are caused by the 
grid that is no longer an optimal representation of the properties 
of the medium. A straightforward solution for these problems is 
updating the grid in such way that the physical properties of the 



medium remain correctly accounted for. By ensuring that the op- 
tical depth between grid points is always close to unity, the grid 
remains optimally suited for the radiative transfer calculation. 

In some cases the removal of grid points according to this 
scheme leads to regions devoid of grid points. An example is 
the photo-ionisation of a cloud of neutral hydrogen gas, which 
causes the optical depth for ionising photons to drop so dramat- 
ically that almost no grid points should be placed if the optical 
depth between grid points has to be of order unity. This extreme 
example leads to different errors in the solution than the ones 
previously described. For example, recombination rates will be 
incorrect as the density in the cloud is no longer resolved by any 
grid point. We circumvent this by imposing a minimum resolu- 
tion below which no grid points are removed. This ensures that 
in every part of the simulation relevant structures are resolved 
and our requirements for accuracy are met. The effect of grid 
point removal and the optimal value for the minimum resolution 
in r ealisti c applications will be further explored in Sects. 15.1.11 
and l533l 

The consequence of preventing undersampling errors is that 
there will remain grid points in the simulation domain between 
which the optical depth is lower than unity. Numerical errors in 
the transport of photons between these grid points are prevented 
by using the transport scheme described in Sect. 12.2.31 As we 
will show in that section, this transport scheme is more expen- 
sive both in computation time as in memory usage compared to 
the other transport schemes. For optimal computation time and 
memory usage it is therefore important to keep the number of su- 
perfluous grid points in regions with low opacity to a minimum. 

2.2. Radiation Transport 

We have shown how the unstructured grid is created on which 
the radiative transfer calculation will be performed. In this sec- 
tion we will show how we can employ the unique properties of 
this grid to efficiently and accurately solve the radiative transfer 
equation. 

During a radiative transfer calculation photons are trans- 
ported from grid point to grid point along the edges of the 
Delaunay triangulation. At every grid point an interaction 
takes place, with the interaction coefficient given by Eq. (fl2b . 
According to the solution of the 1 -dimensional radiative transfer 
equation the number of photons that interacts at this grid point is 

Nact = N in (l - e~ c \ (17) 

where Nj n is the number of incoming photons. The number of 
photons that is not interacting at the grid point is 

N out = N in e- c , (18) 

these photons should continue along their original path. 

In this photon transport scheme there is no difference be- 
tween a grid point that is a source and a grid point that is not. A 
source is simply defined as a grid point that sends photons to all 
its neighbours, which has no influence on the number of compu- 
tations involved. Thus, the SimpleX method does not scale with 
the number of sources. 

After interaction at a grid point there are three ways of pho- 
ton transport, depending on the opacity regime and the physical 
process at hand. 

2.2.1. Scattering processes 

In the case of an isotropic scattering process, or the absorption 
and re-emission of photons at a grid point, the outgoing pho- 
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Fig. 4. Two-dimensional examples of the three modes of transport with which photons are transported from one grid point to another. 
Red arrows indicate incoming photons, blue arrows outgoing photons. Transport in case of a scattering process is shown on the left 
hand side. Incoming radiation loses all memory of the initial direction and is sent to all neighbours of the vertex. Ballistic transport 
is shown in the centre plot. The photons are redistributed to the 2 most straightforward neighbours of the vertex, with respect 
to the Delaunay edge of the incoming photons. On the right hand side direction conserving transport is shown. The photons are 
redistributed to the 2 most straightforward neighbours with respect to the Delaunay edge that is associated with the direction bin the 
incoming photons are in. The outgoing photons stay in the same direction bin and have thus a memory of their original direction. 



tons have no memory of their original direction. In SimpleX 
these photons are isotropically redistributed to all neighbouring 
vertices, as depicted on the left hand side of Fig. |H This type 
of transport does not increase the number of computations in- 
volved, SimpleX is therefore ideally suited to simulate scattering 
processes. 

For the application to ionising radiation it is straightfor- 
ward to include the diffuse radiation from recombining atoms 
in the simulation. Hydrogen ions that capture an electron di- 
rectly to the n = 1 state emit photons capable of ionising 
other atoms. Most radiative transfer methods do not include 
this radiation but instead adopt th e on-the-spot approximation 
(e.g. lOsterbrock & Ferlandl (|2006)), assuming these photo ns to 
be absorbed close to the emitting atom (see also Sect. 14.2b . 
Even though the validity of this approach is not well established 
( Ritzerveld 2005) we will use the on-the-spot approximation for 
all the tests presented in this paper, in order to make a clean 
comparison between our results and the analytic solution or the 
results of other codes that all use the on-the-spot approximation. 

2.2.2. Ballistic transport 

The N out photons in Eq.[l8]that are not interacting at a grid point 
should continue travelling in their original direction. However, 
on the unstructured grid there is no outgoing Delaunay edge in 
the same direction as the incoming edge. This is solved by split- 
ting the photon packets into p equal parts and dividing these 
packets among the p most straightforward neighbours with re- 
spect to their incoming direction. The total number of photons 
is conserved in this process. Tests indicate that if we choose p 
equal to the dimension of the problem d, the solid angle that is 
represented by one Delaunay edge of the emitting vertex is best 
preserved. In other words, a source that sends out photons in all 
directions will fill the entire 'sky' with radiation. However, on 
the unstructured grid it might be the case that one of the most 
straightforward neighbours lies outside a cone of 90 degrees. 
To prevent photons form travelling backwards, we exclude those 



neighbours. In Fig.|4lthe centre plot shows an example of ballis- 
tic transport. 

The advantage of this transport scheme is that the most 
straightforward directions have to be calculated only once, at 
the start of the simulation. As long as the grid doesn't change 
these directions do not have to be recalculated. One important 
disadvantage of this transport scheme is that there is no memory 
of the original direction of the photons. At every interaction the 
outgoing direction was computed from the incoming direction 
in that step, so deflections from the original direction can add 
up, causing numerical diffusion to dominate after approximately 
5 interactions (KPCI09). As long as the mean free path of pho- 
tons is smaller than 5 Delaunay edges, this numerical diffusion 
has no influence on the results, since photons will be interacting 
with the medium before the diffusion becomes dominant. This 
means that during the grid calculation we have to be careful that 
the interaction coefficient c in Eq. (Tf2l) is close to one or larger. 
However, this may lead to too severe constraints on the num- 
ber of grid points that can be placed in optically thin regions. 
Therefore, a different type of transport can be employed in opti- 
cally thin media. 

2.2.3. Direction conserving transport 

If the interaction coefficient c in Eq. (Tf2l) is smaller than one, 
it is no longer sufficient to determine the direction of the pho- 
tons from the direction in the previous step, but a memory of the 
initial direction of the photon is needed. If every photon would 
remember its initial direction and at every interaction point the 
next interaction point would have to be calculated from this di- 
rection, the computation time in optically thin regions would 
grow unacceptably. Instead, the original direction is preserved 
by confining photons to solid angles corresponding to global di- 
rections in space. Unless interacting with the grid, photons stay 
in the same solid angle as they travel along the grid. Even though 
the direction of the photons is now decoupled from the grid, the 
photons still travel along the edges of the triangulation in the 
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Fig. 5. Example of the photon path deviating from a straight line 
due to the unstructured grid. Photons that are travelling in the di- 
rection of the Delaunay edge coming in from the left, should be 
travelling in a straight line along the dotted blue line. However, 
as this is impossible on the unstructured grid, photons travel 
along the Delaunay edges closest to their original direction. This 
path is depicted by the red arrows. The total path depicted by 
the red arrows is longer than the length of the blue line, which is 
corrected for by a global factor. 



same manner as during ballistic transport. Direction conserving 
transport is shown on the right hand side of Fig.ffl 

Since photons still travel along the edges of the triangulation, 
the photon path deviates from an exact straight line in which 
the photons should be travelling, see Fig. Even though the 
direction of the photons is preserved, their paths are longer than 
physically possible. In other words, photons travel slower than 
the light speed on the unstructured grid. We solve this problem 
by applying a global correction factor to the distance between 
grid points, thus ensuring photons travel with the correct speed. 

Introducing these global directions on the unstructured grid 
gives rise to preferential directions, one of the problems the un- 
structured grid was meant to solve. By rotating the solid angles 
over random angles in between photon transport, the preferential 
directions disappear. The drawback of this procedure is that it 
makes direction conserving transport computationally more ex- 
pensive than ballistic transport. For the latter, the most straight- 
forward directions are calculated from the grid and thus have to 
be calculated only once, at the start of the simulation. For direc- 
tion conserving transport it is necessary to recalculate the most 
straightforward directions every time the direction bins rotate. 
Another drawback of direction conserving transport is that the 
photons now have to be stored in n direction bins instead of on 
average 16 neighbours. Typicall y, n = 42 gives converged re- 
sults, but as we will see in Sect. 15.21 this depends on the num- 
ber of optically thin grid points the photons traversed. Thus, 
the memory requirements for direction conserving transport are 
higher. 

2.2.4. Combined transport 

The three modes of transport described above are in general 
applied simultaneously during a simulation. Depending on the 
physical process at hand, photons are transported to all neigh- 
bours (diffuse transport), or to the d most straightforward neigh- 
bours (ballistic or direction conserving transport). In regions 
where the optical depth is higher than or close to one, ballis- 
tic transport is used, while in the optically thin regions direction 
conserving transport is applied. 

Of the three modes of transport, direction conserving trans- 
port is computationally the most expensive (see Sect. 13.31 for a 
comparison between the computation time of ballistic and di- 
rection conserving transport). By applying this scheme only in 



the regions where it is necessary, the computation time is dras- 
tically reduced. As mentioned earlier, numerical diffusion starts 
to dominate in ballistic transport after approximately 5 steps. A 
first guess would therefore be to switch from ballistic to direction 
conserving transport when the optical depth after 5 interactions 
is less than one. That way, we are sure that the majority of pho- 
tons does not take more than 5 steps in ballistic transport. The 
influence of the optical depth at which is switched in a realistic 
simulation is studied in more detail in Sect. 15.1.31 Another way 
to reduce the computation time is by applying the grid dynam- 
ics scheme from Sect. 12.1.41 Removing superfluous grid points 
in the low opacity regime limits the number of vertices at which 
direction conserving transport is performed. 

In combined transport we need to convert from one trans- 
port scheme to another. This is straightforward because every 
Delaunay edge of an optically thin vertex is associated with a 
solid angle in a global direction. When this vertex sends photons 
to an optically thick vertex the photons are transported along the 
Delaunay edges, so the optically thick vertex stores the photons 
according to the Delaunay edge associated with the solid angle. 
In the opposite case, when an optically thick vertex sends pho- 
tons to an optically thin vertex, the photons are converted to the 
solid angles associated with the Delaunay edge along which the 
photons were sent. 



3. Parallellisation strategy 

Even though the radiative transfer scheme presented in the pre- 
vious sections is computationally efficient, in order to do large 
simulations involving a very high number of grid points it is nec- 
essary that the algorithm can run in parallel on distributed mem- 
ory machines. This will not only reduce the computation time 
involved, it also reduces the amount of memory needed at each 
processor to store the physical properties of the grid points. The 
transport algorithm described in Sect. l2.2l has the advantage that 
it is local: the only information needed to do a radiative transfer 
calculation is stored at the neighbours of the vertex. This makes 
the method relatively easy to parallellise. By choosing a smart 
domain decomposition we can minimise the number of commu- 
nications involved. 

3.1. Domain decomposition 

The computation time of a SimpleX calculation is independent 
of the number of sources, it is therefore sufficient to have a do- 
main decomposition that assigns every processor an approxi- 
mately equal number of grid points. Dividing space into equal 
volumes and assigning each volume to a processor is not suffi- 
cient because the number of points in each volume may differ 
dramatically due to the adaptive grid. We therefore use a do- 
main decomposition based on the space-filling Hilbert curve, 
whic h is also employed in other methods without a re gular 
grid (IShirokov & Bertschingedl20Q5l:ISpringelll2QQ5Ll2QlQl) . The 
Hilbert curve is a fractal that completely fills a cubic rectangular 
volume. A Hilbert curve is uniquely defined by its order m and 
its dimensionality d, filling every cell of a ^-dimensional cube 
of length 2 m . The following properties of the Hilbert curve are 
beneficial when using it for domain decomposition: 

- Locality: points that are close along the ID Hilbert curve are 
in general also close in 3D space. 

- Compactness: a set of cells defined by a continuous section 
of the Hilbert curve has a small surface to volume ratio. 
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- Self-similarity: the Hilbert curve can be extended to arbitrar- 
ily large size. 

The first two properties minimise the number of communications 
between processors, while the third property ensures that we can 
use an arbitrarily large number of cells to determine the domain 
decomposition. 

The first step in the domain decomposition is dividing the 
domain into 2 md equal, regular cells, where d is again the dimen- 
sion and m the order of the Hilbert curve. We then step through 
the cells along the Hilbert curve, counting the number of grid 
points inside each cell until the number of grid points equals the 
total number of grid points divided by the number of processors. 
In this way, every processor approximately holds an equal num- 
ber of grid points, thus dividing the work load evenly, while the 
necessary communications between processors are minimal due 
to the locality and compactness property of the Hilbert curve. 

3.2. Parallel radiative transfer 

The QHull algorithm that is applied to calculate the triangula- 
tion works only in serial. In the case of parallel execution of 
SimpleX, every processor calculates the triangulation of the ver- 
tices that belong to that processor and vertices in a border around 
it belonging to neighbouring processes. This border is used to 
connect the triangulation between different processes. We ensure 
that the border is large enough by using the empty circumsphere 
principle. For all simplices that contain at least one vertex inside 
the domain of the current process, we ensure that the circum- 
sphere of the simplex lies entirely within the boundary around 
the domain. Thus, we are certain that no vertex exists on another 
process that lies inside the circumsphere of this simplex and the 
triangulation is valid. After the triangulation algorithm has been 
applied, every process keeps only the vertices assigned to the 
process and a local copy of vertices assigned to other processes 
that are neighbour to a vertex on this process. These local copies 
are strictly used to send photons from one process to another, no 
physical interaction is taking place. 

The transport scheme described in Sect. l2.2l is local, because 
photons are only transported from one grid point to another. We 
can therefore do a full radiative transfer time step without any 
communication between processes. During a time step, photons 
might be sent to a neighbour of a vertex that does not exist on 
the current process. After each time step, these photons are then 
communicated to the appropriate processes and the cycle starts 
anew. The local copies of vertices are only used for the trans- 
port of photons, all the physical interactions are taking place on 
the process that the vertex is assigned to. Hence, we are certain 
that physical interactions take place exactly once for every ver- 
tex during a radiative transfer time step. 

3.3. Scaling tests 

To check how well the described parallellisation strategy works, 
we have conducted several scaling tests. The test consists of a 
simulation of a single source in a homogeneous medium, simi- 
lar to the set-up of the test in Sect. 15.11 with the difference that 
the simulation is stopped at 10 Myr. No grid points were re- 
moved during these tests. All the simulations were conducted 
using AMD Opteron 246 64Bit CPUs of 2.6 Ghz with 4 GB of 
memory per node. Unfortunately, we had only 8 nodes available 
for these tests, but it will give a general idea of the scaling. Note 
that even though we have used only one source for these scaling 
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Fig. 6. Simulation time as a function of the number of grid 
points. Shown are the computation time of the triangulation 
(black solid curve), ballistic radiation transport (red dotted), 
combined radiation transport (cyan long dashed), direction con- 
serving transport (green short dashed), total simulation time with 
ballistic transport (blue dot- short dashed), total simulation time 
with combined transport (brown short dash-long dashed) and 
total simulation time with direction conserving transport (vio- 
let dot-long dashed). A doubling of the number of grid points 
N means a doubling of the computation time for most parts 
of the simulation. An exception is the triangulation algorithm 
that scales O (N log N) and the combined transport scheme. The 
computation time of the latter depends strongly on the number 
of optically thin grid points and thus on the physics of the sim- 
ulation at hand, but will always be in between the computation 
time of ballistic transport and direction conserving transport. 



tests, the inclusion of more sources would not have influenced 
the timings presented here. 

As a first test we used only one node with an increasing num- 
ber of grid points, shown in Fig. [6] This figure shows that an in- 
crease of the number of grid points N increases the computation 
time linearly for most components of the simulation. Exceptions 
are the triangulation algorithm, which scales O(NlogN) and 
the combined transport scheme. The computation time of the 
combined scheme will always be between the computation time 
of ballistic transport and that of direction conserving transport. 
It depends highly on the number of optically thin grid points. 
For the low resolution simulations, the number of optically thin 
grid points is relatively high and therefore the computation time 
is comparable to the computation time of direction conserving 
transport. Increasing the number of grid points decreases the rel- 
ative number of optically thin grid points and therefore the com- 
putation time of combined transport comes closer to that of bal- 
listic transport in the case of more grid points. Note that this is a 
feature of the set-up that we chose for the scaling test. If a larger 
region of the computational volume would be ionised, the com- 
putation time would be longer. On the other hand, if the ionised 
region were smaller, the computation would be shorter. Fig. [6] 
also shows that the computation time of the entire simulation is 
dominated by the radiation transport and not by the construction 
of the triangulation, since the computation time of the triangula- 
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Fig. 7. Simulation time as a function of the number of processors 
for a constant number of grid points. The number of grid points 
is 128 3 = 2097152. Shown are the computation time of the 
triangulation (slate dashed curve), ballistic radiation transport 
(red dotted), combined radiation transport (cyan long dashed), 
direction conserving transport (green short dashed), total sim- 
ulation time with ballistic transport (blue dot- short dahed), to- 
tal simulation time with combined transport (brown short dash- 
long dashed) and total simulation time with direction conserving 
transport (violet dot-long dashed). Most components of the sim- 
ulation scale very close to linear as the number of processors 
increases. An exception is the triangulation algorithm, due to the 
fact that every processor needs to triangulate extra points that are 
in the boundary between processors. 

tion is one order of magnitude shorter than that of the radiation 
transport. 

Fig. [7] shows the strong scaling properties of the SimpleX al- 
gorithm. We simulate the same physical problem as the previous 
test, but this time the number of grid points is held constant at 
128 3 = 2097152. By increasing the number of processors we 
can analyse the extra work that needs to be done when more 
processors are employed. In the ideal case no extra work would 
have to be done at all, this would result in the black solid curve 
shown in the figure. The only component of the simulation that 
does not follow the linear scaling very well is the triangulation, 
which can be easily understood by the way the triangulation is 
constructed on multiple processors. Every processor has to cal- 
culate the triangulation of the grid points assigned to the proces- 
sor and an additional number of grid points in a boundary around 
this domain. Increasing the number of processors thus means ef- 
fectively increasing the number of grid points that needs to be 
triangulated and therefore the scaling is not as favourable as one 
might hope. However, the computation time of the triangulation 
remains an order of magnitude smaller than the radiative trans- 
fer components that do scale almost linearly, so this presents no 
serious issue. 

Finally, Fig.[8]shows the weak scaling of the algorithm. With 
the same problem set-up as before, we now increase the num- 
ber of processors while keeping the number of grid points per 
processor constant. Ideally, the amount of work per processor 
would stay the same and the computation time would remain 



Fig. 8. Simulation time as a function of the number of processors 
for a constant number of grid points per processor. The number 
of grid points at each processor is 64 3 = 262144. Shown are the 
computation time of the triangulation (black solid curve), bal- 
listic radiation transport (red dotted), combined radiation trans- 
port (cyan long dashed), direction conserving transport (green 
short dashed), total simulation time with ballistic transport (blue 
dot- short dahed), total simulation time with combined trans- 
port (brown short dash-long dashed) and total simulation time 
with direction conserving transport (violet dot-long dashed). The 
computation time of the radiation transport components shows 
a marginal increase as the problem size gets bigger due an in- 
crease in the number of communications involved. The compu- 
tation time of the triangulation increases more for the same rea- 
son as before, that more boundary points have to be triangulated 
when the number of processors increases. 



constant. Within a few percent the computation time for the ra- 
diation transport components remains constant, it only increases 
marginally due to an increase in the amount of communication 
necessary. As with the strong scaling, the triangulation algorithm 
needs to triangulate more boundary points as the number of pro- 
cessors increases and therefore the computation time increases 
with number of processors. However, the number of boundary 
points is small compared to the total number of grid points that 
needs to be triangulated, therefore the curve of the computation 
time flattens with an increase of the number of processors. 



4. Radiative transfer of ionising radiation 

Having laid down the basics of the transport mechanism of 
SimpleX, in this section we will describe the application to ion- 
ising radiation. Generally, the input for the radiative transfer cal- 
culation will be given by a hydrodynamics simulation that can 
be in the form of a regular, AMR or SPH grid, from which the 
Simp leX grid is calculated according to the recipe described in 
Sect. 12.11 Photons are then transported on this grid as described 
in Sect. 12.21 Our aim is to apply SimpleX to cosmological ap- 
plications, we therefore start this section with the cosmological 
radiative transfer equation. However, the method can be used for 
all applications that require the transport of ionising photons. 
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4.1. Cosmological radiative transfer equation 

The cosmological radia t ive transfer equation reads (e.g. 
iGnedin & Ostrikerl (ll997l) : lNorman et alJ (1 1998b ) 



Idly fi VIy H{t) I dly . . 

+ Z V- 3Iy = J V - (Xyly, 

c ot a c \ av 



(19) 



where I v = I(h, x, t, v) is the specific intensity with frequency v 
along the unit propagation direction vector h at position x and 
time t, j v (x, h) is the emission coefficient, and a v (x, h) is the 
extinction coefficient. The latter includes the scattering coeffi- 
cient a s v cat (x, h) and the pure absorption coefficient a^ bs (x, h). 
Furthermore, H(t) = a I a is the time-dependent Hubble parame- 



ter, and a 



l+Z 



is the ratio of cosmic scale factors between 



photon emission and present time t. If the mean free path of 
photons is much smaller than the horizon, the classical transfer 
equation is a valid approximation: 



1 dl 

-— ^ + h • VI V - j v (x, h) - a v ix, h)I v . 
c ot 



(20) 



This approximation holds fairly well during the beginning of 
reionisation. However, care must be taken when ionised bubbles 
start to overlap and the photon mean free path increases dramat- 
ically. In this case the expansion of the Universe does have to be 
taken into account. 

One more simplification can be made if j v and a v change on 
time scales larger than the light crossing time in the simulation 
volume. In this case the time-dependence can be dropped: 



Q • Vly = j v (x, h) - a v ix, h)Iy. 



(21) 



It is this equation that SimpleX solves. In most astrophysical 
fluid flows this approximation holds fairly well as long as time 
scales over which the radiative transfer equation is solved are 
sufficiently short. However, one must take care that this equa- 
tion implicitly assumes an infinite speed of light, which might 
re sult in unphysical behaviour. For example, as was pointed out 
bv lAbel et all (1 1999) close to ionising sources an ionisation front 
might travel faster than the speed of light. In the tests presented 
in this paper, we have checked the validity of this approximation. 



4.2. Ionisation processes 

Photons with frequencies above the Lyman limit vo can be ab- 
sorbed by neutral hydrogen atoms. The n umber of photoionisa- 
tions per unit time per hydrogen atom is (lOsterbrock & Fer land 
2006): 



H Jvo hy 



cr H (v)dv, 



(22) 



where J v is the mean intensity and <th(v) is the ionisation cross 
section for H by photons with energy hv > Uvq. This cross sec- 
tion can be approximated by 



crn(v) 



V v > vo, 



(23) 



with reference value Ao = cr H (y ) = 6.3 • 10" 18 cm 2 . The in- 
verse process is the recombination of electrons with the ions. 
The number of recombinations per hydrogen atom per unit time 
is: 



where au(T) is the recombination coefficient of hydrogen atoms 
and n e the electron density. The evolution of the number density 
of ionised atoms at a certain point in space can then be written 

as 



dt 



(25) 



where nuih ^hi and nu are respectively the number density of 
ionised hydrogen atoms, neutral hydrogen atoms and the total 
number of atoms. In this equation we have neglected collisional 
ionisations, but they can be included trivially by adding an extra 
term to the equation. Relevant time scales for these processes are 
the photoionisation time scale ti on = 1 /T and the recombination 
time scale t rec = l/R. 

When a hydrogen ion captures an electron directly to the 
ground level, radiation is emitted with a frequency above the 
Lyman limit. In almost all radiative transfer codes, it is assumed 
that this radiation is absorbed close to the emitting atom, the 
so-called on-the-spot approximation. One can therefore ignore 
recombinations to the ground level, as they are cancelled out by 
the emitted radiation, and u se the correspondi ng recombination 
coefficient (Xb(T). However. iRitzerveldl (120051) showed that, de- 
pending on the density distribution, if the source produces ra- 
diation just above the Lyman limit, this approximation is not 
valid. Most radiative transfer codes are incapable of including 
this diffuse recombination radiation since effectively every grid 
cell that is ionised becomes a source. In the SimpleX algorithm it 
is straightforward to include the diffuse recombination radiation 
self-consistently. However, the analytic solutions and the results 
from other codes with which we will be comparing the SimpleX 
results in the next section all rely on the on-the-spot approxima- 
tion. Therefore, we will use the on-the-spot approximation for 
all tests presented in this paper. 

4.3. Assigning sources 

On the SimpleX grid, a source is defined as a grid point that sends 
photons to all of its neighbours. For all the tests presented in this 
paper, we use a single frequency for the photons. The luminosity 
of the source is obtained by integrating the source spectrum S v 
divided by the energy of each photon over frequencies higher 
than the Lyman limit vq: 



= r—dv 



(26) 



By integrating the source spectrum over frequency, we neglect 
the influence of the spectrum on the interaction of the photons 
with the medium. We compensate for this b y using a mean opac- 
ity representation fsee lMihalas & Mihalai (11984 ): 



(T = L 



-1 



Jo hv 



(27) 



The number of photons that a source sends out at every step is 
determined by the source luminosity and the radiative transfer 
time step At rt . 

4.4. Interaction 

Photons that were sent out by a source travel in one radiative 
transfer time step from grid point to grid point, where an interac- 
tion with the medium takes place. The photons travel a distance 
AL between grid points, and thus encounter an optical depth 



Ru = n e an(T), 



(24) At = n m o-AL = (1 - x)ft H AL, 



(28) 
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in which x is the ionised fraction inside the Voronoi cell through 
which the photon travels. Except for the ionised fraction, all 
these quantities have been calculated during the creation of the 
grid. Eq. (|28t is equivalent to the interaction coefficient in Eq. 
JT2l) . Thus, if the incoming number of ionising photons is Nj n , 
the number of ionising photons that is absorbed at this grid point 
is 



Nabs = N in (l - e~ AT ), 



(29) 



and the number of ionising photons that is propagating onwards 
is 



Nam = N in e 



(30) 



The number of photons that are absorbed ionise the medium, 
thereby changing the local ionised fraction. As the medium gets 
ionised, the optical depth at this grid point changes, which means 
we should either use the direction preserving scheme at this grid 
point, or remove grid points to preserve the relation between op- 
tical depth and Delaunay line length, as described in Sect. 12.21 
and Sect. [2X1 

4.5. Solving the photoionisation rate equation 

Having established the number of photons that is available for 
ionising the cell, it is straightforward to convert this to a pho- 
toionisation rate and solve Eq. (l25t . However, care must be taken 
since by doing this we implicitly assume that the neutral den- 
sity stays constant during a radiative transfer time step. This is 
only true for At n «: ft on and At n <^c We t herefore adopt the 
scheme described in iPawlik & Schayd (120081) and subcycle the 
rate equation on time steps A£ C h em < Af rt assuming that the ion- 
ising flux is constant during a radiative transfer time step. This 
ensures photon conservation even if the radiative transfer time 
steps are large. The rate equation at time t C h Qm e (t n , t n + Af rt ) is 
then 

^^(^chem) „(^chem)"T , (/chem) A + T/ ,(/chem) ...fchem) /T\ 1 \ 

an mi ~ n Ul A H A ^chem-^e n mi a ^ ) A ^hem, pi) 

where the photoionisation rate at t c ^ m is given by 



r(^chem) T" 
H " 1 H 



1— T^chem^ \ 



1 



(/chem ) 



HI 



where Th and r are the photoionisation rate and optical depth 
at the beginning of the subcycling and r (?chem) = t«||j /aihi- 
By denning the photionisation rate in this way, the ionising 
flux in the cell is constant during the radiative transfer time 
step. This subcycling scheme becomes computationally expen- 
sive when A^chem <^ Af rt , but photoionisation equilibrium is gen- 
erally reached after a few subcycles. It is then no longer neces- 
sary to explicitly integrate the rate equation, but instead use the 
values of the preceding subcycle step. This way of subcycling en- 
sures photon conservation even for large radiative transfer time 
steps. 



4.6. Time stepping 

As we discussed in Sect. 14.11 we are primarily interested in 
solving the time-independent radiative transfer equation, which 
means that the speed of light is assumed to be infinite. In all 
the tests presented in this paper, photons are moved only one 
Delaunay edge during a radiative transfer time step, thus inter- 
acting only at one grid point in that time step. We therefore have 



to be careful that the radiative transfer time step At n is suffi- 
ciently small to satisfy the time-independent transfer equation. 
In the limit that At n goes to zero, the time it takes for photons 
to leave the simulation box goes to zero as well, satisfying the 
condition of infinite speed of light. For all tests presented in this 
paper, we have checked that the radiative time step is sufficiently 
small to be in agreement with this limit. 

The subcycling scheme that is used to calculate the evolution 
of the neutral fraction at a grid point during a radiative transfer 
time step allows for much larger time steps than needed to satisfy 
the time-independent transfer equation. This is very useful in 
simulations where the photons are allowed to travel more than 
one Delaunay edge per time step, for example in case one needs 
to solve the time-dependent transfer equation. However, this was 
not done for the tests presented in this paper. 

5. Tests 

In order to test the accuracy of the new SimpleX algorithm, 
we have performed several te sts that were par t of the Radiative 
Transfer Comparison Project (llliev et al .1120061) . The original im- 
plementation of SimpleX was part of this project but only did 
some of the tests. This gives us the opportunity to show the dif- 
ferences between the two versions in these tests, and show the 
behaviour of the new algorithm in a shadowing test that wasn't 
originally done by SimpleX for the Comparison Project. For 
comparison purposes, in all tests presented here we adopt the 
on-the-spot approximation. 

The increasing realism of the tests presented here allows us 
to highlight the different improvements of the method. The first 
test shows the importance of the direction conserving transport 
scheme to accurately account for the ionisations in the regions 
with low opacity and quantifies the optical depth at which it is 
necessary to switch from ballistic to direction conserving trans- 
port. The second test enables us to quantify the number of di- 
rection bins necessary to account for shadowing behind a dense 
cloud in direction conserving transport. Finally, in the third test 
we can assess the importance of the new sampling routine and 
study the minimum resolution required when dynamic grid up- 
dates are applied. 



(32) 5. 7. Test 1: Isothermal Hn region expansion 



One of the few problems in radiative transfer that has a known 
analytical solution is the Hn region expansion in a homoge- 
neous medium. A steady monochromatic source emits N y pho- 
tons per second of frequency hv = 13.6 eV into an initially 
neutral medium with constant gas density hr- In equilibrium, 
the number of photons emitted by the source is balanced by the 
number of photons absorbed due to recombinations in a spher- 
ical volume. The radius at which equilibrium is reached is the 
Stromgren radius, given by 



3Nv 



1/3 



rs = lw^J ■ (33) 

Assuming an infinitely thin ionisation front and a fully ionised 
inner region, the ionisation front radius and velocity as function 
of time are 



n = r s 
and 



1/3 



v 7 = 



3 *rec (1 -e-^rec) 2/3 ' 



(34) 



(35) 
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shown as the black solid lines in Fig. [9j 

We can improve on this by dropping the assumption of a 
fully ionised inner region of the Stromgren sphere and calculat- 
ing the neut ral and ionised fraction as a function of radius by 
solving (e.g. lOsterbrock & Ferlandl (120061) ) 



Anr 2 J 



dvN y (v)e-^(Tn(y) = J(r)n 2 R a B (T). 



(36) 



By using the commonly employed definition of the position of 
the ionisation front as the radius at which x = 0.5, solving this 
equation by direct integration gives us a second way of obtaining 
the Stromgren radius. This yields a slightly different ionisation 
front position than obtained in Eq. (l34l) . We show the solution 
obtained from directly integrating Eq. (l36t as the black dashed 
line in Fig.0 

The analytical solutions thus obtained make this test ide- 
ally suited to test the different transport types of SimpleX de- 
scribed in Sect. 12.21 We therefore performed this test problem 
with ballistic, direction conserving and combined transport to in- 
vestigate the behaviour of the specific schemes. The parameters 
for this test are as follows. The computational box has length 
L = 13.2 kpc, the gas number density is = 10~ 3 cm~ 3 , the 
temperature of the gas is T = 10 4 K. A source is placed in the 
centre of the box, emitting N y = 5 - 10 48 ionising photons s" 1 . 



For these parameters, t vt 



3.86 • KPs = 122.4 Myr and 



rs = 5.4 kpc. The total simulation time is 500 Myr » At r e o. Not e 
that this test differs slightly from Test 1 in llliev et al.1 00061) . 



where the computational volume is smaller and the source is 
located in the corner of the computational box. Except where 
noted, a resolution of 64 3 grid points and a time step of 0.05 Myr 
is used for this test. The grid on which we will perform this test 
is constructed by using the recipe described in Sect. 12.11 by us- 
ing a homogeneous Poisson process to place the grid points. This 
introduces more shot noise compared to using a glass-like distri- 
bution, in which the point process is modified to make the points 
avoid one another. However, this is the same procedure we will 
apply for inhomogeneous density distributions using Eq. (fT3t . 
so in order to get a good understanding of the limitations of the 
method, we choose to use the Poisson process over a glass-like 
distribution for this test. 



5.1.1. Ballistic transport 

The single mode of transport in the original SimpleX algorithm 
was ballistic transport, described in Sect. I2.2.2I This mode of 
transport was designed for regions where r > 1 between grid 
points. However, as the medium gets ionised during the simu- 
lation, the optical depth between grid points becomes so small 
that it is no longer correct to transport photons in this way. As 
described in KPCI09, a photon transported ballistically loses all 
memory of its initial direction after approximately 5 steps on the 
grid. Therefore, using ballistic transport in the highly ionised 
inner region of the Stromgren sphere introduces numerical dif- 
fusion. 

The numerical diffusion does not influence the position of 
the ionisation front severely. As is shown in Fig. the red dot- 
ted line representing ballistic transport follows the numerical so- 
lution (the black dashed line) very closely, the error at the end 
of the simulation time is approximately 1 percent. The inner 
structure of the ionised region will be wrong, however, as we 
expect the numerical diffusion to dominate in the inner region 
of the Stromgren sphere, were a large number of steps in an 
optically thin region needs to be taken, instead of close to the 
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Fig. 9. The position, relative error and velocity of the ionisation 
front for Test 1 . The black solid curves represent the analytic so- 
lutions of Eq. (l34l) and Eq. (|35K while the black dashed curve 
represents the results of directly integrating Eq. (f36b . Shown in 
colour are the results of SimpleX simulations with ballistic trans- 
port only, where the red curves represent a simulation with a 
static grid and the violet, blue and green curves simulations with 
a dynamic grid with minimum resolutions of 48, 32 and 16, re- 
spectively. The position of the ionisation front is within 1 % of 
the analytical solution, although the effects of the limited reso- 
lution are clearly visible in the runs with a low minimum resolu- 
tion. 



ionisation front. In Fig.[l0l the spherically averaged neutral and 
ionised fractions are plotted as a function of distance from the 
source. The numerical diffusion in the inner region results in a 
too low neutral fraction. The reason for this is that the source 
photons quickly become diffuse and therefore, instead of travel- 
ling straight to the ionisation front, stay longer in the inner parts, 
thus cancelling more recombinations and causing a lower neutral 
fraction than expected from the analytical solution. 

One possible solution to this problem is remov ing grid points 
that have too low optical depths (see Sect. 12.1.41) . In Fig. [9] and 
Fig. [TO] the results are shown for simulations where grid points 
are removed until a certain minimum resolution is reached. Fig. 
[TO] shows that the effect of numerical diffusion on the inner struc- 
ture of the ionised region is lessened by the removal of grid 
points. The neutral fraction comes closer to the analytic solu- 
tion as the minimum resolution decreases and the equilibrium 
position of the ionisation front becomes slightly more accurate. 
However, the slightly more accurate equilibrium results come at 
a cost. In Fig. [9] we see that the lower resolution of 16 and 32 
in the inner region causes the ionisation front position to devi- 
ate more than 5 percent from the analytic solution at early times, 
even though the equilibrium solution is accurate. Also, the spher- 
ically averaged equilibrium neutral fraction shows some severe 
artefacts due to the low resolution, most pronounced close to the 
source. Therefore, we conclude that only removing grid points 
with low optical depth is not a viable remedy against numerical 
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Fig. 10. Spherically averaged neutral and ionised fractions as 
function of the radial distance from the source after 500 Myr for 
Test 1 . The black solid curve represents the result of directly in- 
tegrating Eq. d36b . Shown in colour are the results of SimpleX 
simulations with ballistic transport only, where the red curve 
represents a simulation with a static grid and the violet, blue 
and green curves represent simulations with a dynamic grid with 
minimum resolutions of 48, 32 and 16, respectively. Ballistic 
transport in the ionised regions results in numerical diffusion, 
therefore the neutral fraction is too low in the inner region, as can 
be seen from the dotted red curve. Removing grid points during 
the simulation, so essentially keeping the optical depth between 
grid points constant during the simulation, alleviates the diffu- 
sion problem, but introduces other errors as a result of the low 
resolution. The minimum resolution imposed gives a handle on 
how to control the errors stemming from numerical diffusion and 
a too low resolution. 

diffusion, since the low resolution needed in the ionised regions 
causes severe noise in the equilibrium solution. 

5.1 .2. Direction conserving transport 

The numerical diffusion in ballistic transport is caused by the 
loss of direction of t he pho tons after a number of interactions at 
grid points. In Sect. I2.2.31 we described how we can cure this 
problem by defining solid angles in which the photons travel. 
The number of solid angles is a measure for the accuracy of 
the direction conservation. In Fig. [TT] the results of direction 
conserving transport with 21, 42, 63 and 84 direction bins are 
shown. From this plot we can see that the neutral fraction in the 
inner part of the Stromgren sphere follows the analytical solution 
accurately if we use 42 direction bins or more. We can therefore 
conclude that the direction conserving transport scheme is an 
excellent solution for transporting photons in the optically thin 
regime. Since the difference between 42 direction bins and 63 
and 84 bins is negligible, we are justified in using 42 direction 
bins when direction conserving transport is used in this test. In 
the next tests we will show the influence of the number of direc- 
tion bins on the shadowing properties of the algorithm. 

Fig. [TT] also shows that with the use of fewer than 42 di- 
rection bins a small error from the numerical diffusion remains 
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Fig. 11. Same layout as Fig. [TUJ Shown in colour are the re- 
sults of SimpleX simulations with direction conserving transport, 
where the number of directions in which the photons are stored 
is 21, 42, 63 and 84. Clearly, the numerical diffusion of which 
ballistic transport suffers is absent if 42 directions or more are 
used to store the photons. 

present in the simulation, the red dotted line representing the 
simulation using 21 direction bins deviates slightly from the an- 
alytic solution. The angular sampling does improve with the use 
of 21 direction bins compared to the average number of neigh- 
bours of a vertex (approximately equal to 15.54 in this point dis- 
tribution), but the number of bins is too small to prevent photons 
from deviating from their original path. 

The accurate transport of photons in the optically thin regime 
comes at a price. The direction conserving transport is compu- 
tationally more expensive than ballistic transport, due to the fact 
that the direction bins need to be associated with every outgoing 
Delaunay line along which the photons are transported. To pre- 
vent preferential directions on the grid, the direction bins need 
to be rotated randomly after every time step, which causes ad- 
ditional computational overhead. Even though the transport of 
photons itself is almost as fast as with ballistic transport, the 
calculation of the grid properties takes more time. This extra 
computational cost can be reduced by combining both transport 
modes. 

5.1.3. Combined transport 

The introduction of direction bins prevents numerical diffusion 
in the optically thin regime but adds some computational over- 
head. Because the numerical diffusion is only present in the opti- 
cally thin regime, we can speed up the calculation by combining 
ballistic and direction conserving transport in such a way that at 
high and moderate optical depth the faster ballistic transport is 
used, while vertices in the low optical regime employ the direc- 
tion conserving mode of transport. This results in a computation 
time that is significantly faster than direction conserving trans- 
port in typical cosmological applications. 

The optical depth at which is switched from one mode of 
transport to the other is an important parameter. If it is too low, 
ballistic transport is done in regions with too low an optical 
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Fig. 12. Same as Fig. [TO] The coloured graphs represent the re- 
sults of combined transport with different optical depts at which 
is switched from ballistic to direction conserving. If the switch 
is set too low, at r = 0.01 and r = 0.1, the numerical diffusion is 
not completely absent from the simulation. If we choose r > 0.5, 
the results are the same as with completely direction conserving 
transport. 



depth, causing numerical diffusion. If it is too high, direction 
conserving transport is done in optically thick regimes, causing 
unnecessary computational overhead. In Fig. [T2j the influence 
of the optical depth at which is switched is shown. If the conver- 
sion from ballistic to direction conserving transport happens at 
r = 0.01 and r = 0.1, we can see in this plot that there is still 
numerical diffusion in the inner region of the Stromgren sphere. 
However, if the switch is made at r > 0.5, the difference between 
fully direction conserving and combined transport disappears. 
To be on the safe side, we use in the remaining tests a switch at 
r = 1.0. This way, we are sure that the numerical diffusion is 
completely absent in the simulations. 

In order to get a good understanding of the behaviour of the 
new SimpleX algorithm, we used the same time step of 0.05 Myr 
and the same resolution of 64 3 grid points in all previous simu- 
lations. Using the combined transport scheme with the fiducial 
value of r = 1.0 for the switch between combined and direction 
conserving transport, we can proceed to find how large the in- 
fluence of both the time step and the resolution is. In Fig. [13] the 
position of the ionisation front as a function of time is plotted for 
different simulation time steps At. Since photons can travel only 
from one grid point to another during a time step, large time 
steps show unphysical behaviour, especially when the equilib- 
rium solution has not been reached yet. From Fig. [13] we see that 
at a time step of At = 1 Myr the ionisation front is behind the 
analytical solution at early times, but it gives the correct equi- 
librium solution. Time steps smaller than that retrieve the ana- 
lytical solution to within 1%. A time step of 10 Myr gives an 
ionisation front position that is behind even at equilibrium, be- 
cause for these large time steps photons are travelling slower 
than the speed of light. Another reason is that in order to pre- 
vent preferential directions on the grid, the direction bins need 
to be randomly rotated every time step. If the time step is too 
large, there are not enough rotations to avoid these preferential 
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Fig. 13. Same as Fig. [9] The coloured graphs represent the dif- 
ferent time steps of At = 0.01, 0.1, 1 and 10 Myr. 



directions, resulting in an incorrect equilibrium solution. This 
issue could be solved by letting photons travel more than one 
Delaunay edge in a time step. However, in that case the only dif- 
ference between the simulations would be the time at which the 
physics is be evaluated, which is not what we were interested in 
for this comparison. 

Using a time step of 0.05 Myr, Fig. [14] shows the effect of 
the resolution on the neutral and ionised fraction as function of 
radius and compares this to the results of other codes in the 
Radiative Transfer Comparison Project, shown as the shaded 
area. For all simulations, the ionisation front is at the same posi- 
tion. For low resolution, the neutral fraction in the inner part of 
the Stromgren sphere is noisy, but still lies within the range of 
solutions reported by other codes. Higher resolution simulations 
are less noisy and converge to the same solution. 

5.2. Test 2: Shadowing behind a dense cloud 

The second test in the Radiative Transfer Comparison Project 
examines ionisation front trapping in a dense clump and the for- 
mation of a shadow. Unfortunately, this test involves a plane- 
parallel wave-front, which is something that is difficult to impose 
in SimpleX. The reason for this is that a plane-parallel wave as 
described in this test represents a preferential direction, which is 
exactly what we try to avoid in our method. It is possible to alter 
the method to produce a plane-parallel wave. However, it would 
then be unclear how much this alteration affects the shadowing 
properties in other applications where this alteration is not used. 

We therefore chos e to do a different shad owing test, similar 
to the one described in Me llema et al.l (120061) . This test involves 
the irradiation of a dense clump by a point source. The test set- 
up is similar to the test described in the previous section, except 
that we put a dense, uniform, rectangular slab with size 855 by 
2138 by 2138 kpc at a distance of 1096 kpc from the source. The 
density contrast between the homogenous environment and the 
clump is ftciump/ftout = 200. The ionisation front will be trapped 
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Fig. 14. Same as Fig. [TO] The coloured graphs represent simu- 
lations with a different resolution of 16 3 , 32 3 , 64 3 and 128 3 . The 
shaded area represents the range of neutral and ionised fractions 
found by the cod es in the Radiative Transfer Comparison Project 
(llliev et al. 2006). All SimpleX simulations lie within the shaded 
area, showing that even at very low resolution SimpleX gives ac- 
ceptable results. 



inside the dense clump and a shadow should form behind the 
clump. 

This test provides an excellent means to study the number 
of direction bins in which the photons have to be stored at op- 
tically thin vertices to ensure proper shadowing. In Fig. [15] the 
results of this test at a resolution of 64 3 grid points is shown. At 
this resolution the effect of different angular samplings is neg- 
ligible, the shadow that is cast has the same width for both the 
simulation involving 42 direction bins and 84 direction bins. The 
ionisation front penetrates approximately 5 cells into the region 
shadowed by the clump. This is a result of the fact that pho- 
tons are transported along the d most straightforward Delaunay 
edges of the triangulation. Some of these edges are pointing into 
the region shadowed by the clump, thereby ionising the medium 
inside the shadowed region. By restricting the angle in which 
to look for the d most straightforward neighbours (which is 90° 
in these simulations, see Sect. 12. 2K the sharpness of the shadow 
will gradually increase. However, this leads to an ionised region 
that is no longer spherical in the regions that are not shadowed, 
because the source no longer 'fills' the entire sky with radiation. 

Fig. [TO] shows the same test at a resolution of 128 3 grid 
points. At this resolution the effect of using 84 direction bins 
instead of 42 is profound. While with 84 direction bins photons 
penetrate approximately 8 cells into the shadowed region, with 
42 direction bins this is approximately 15 cells. The shadow cast 
in case of 84 direction bins is comparable to the shadow in the 
simulation using 64 3 cells, the width being governed by the d 
most straightforward directions. This is not the case when 42 
direction bins are used to store the photons at optically thin ver- 
tices. There is a small amount of diffusion associated with every 
rotation of the direction bins. At a resolution of 128 3 grid points, 
the photons undergo so many rotations on their path that this 
diffusion starts to be visible in the shadowing properties of the 
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Fig. 15. Test 2 results for a resolution of 64 3 grid points. Shown 
are slices along the xy-plane, xz-plane and yz-plane, as indi- 
cated. The dense clump is shown in black. The contours rep- 
resent the position of the ionisation front at the end of the simu- 
lation, t = 500 Myr. The grey contours represent the simulation 
where the photons are stored in 42 direction bins at optically 
thin vertices, while the black contours represent the simulation 
where 84 direction bins were used. At this resolution, the in- 
creased angular resolution of the photons in the simulation with 
84 direction bins does not have a significant impact on the shad- 
owing. The ionisation front penetrates approximately 5 cells into 
the region shadowed by the clump. The slice through the x-axis 
shows that the ionisation front outside the shadow is not affected 
by the dense clump. 

method. This unwanted numerical scattering is cured by using 
more direction bins. 

This shadowing test has shown that the direction conserving 
transport scheme at optically thin vertices provides the means 
to cast shadows behind dense clumps. As a result of the choice 
to send photons along the edges of the Delaunay triangulation, 
SimpleX will never be able to reproduce the infinitely sharp shad- 
ows that ray tracing methods produce. The random rotation of 
the direction bins at every radiative transfer time step that is nec- 
essary to avoid preferential directions on the grid causes a small 
amount of numerical diffusion that shows up when photons tra- 
verse many optically thin vertices. We have shown that by in- 
creasing the number of direction bins this problem is cured. This 
test shows the importance of removing redundant optically thin 
vertices in highly ionised regions, as described in Sect. 12.1.41 to 
minimise the number of optically thin vertices that need to be 
traversed. 

5.3. Test 3: Multiple sources in a cosmological density field 

The final test we conducted is closest to our intended application, 
that of a cosmological density field. It is for this kind of problem 
that the SimpleX method has been primarily designed, since the 
test consists of multiple sources in a density field with a large 
dynamic range. The initial conditions are given by a time slice at 
z = 9 from a cosmological N-body and g as-dynamic simu lation 
using the cosmological PM+TVD code dRvu et al.lll993h . The 
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Fig. 16. Same as Fig.[l5]but with a resolution of 128 3 grid points. 
At this resolution, the increased angular resolution of the pho- 
tons in the simulation with 84 direction bins does have a sig- 
nificant impact on the shadowing. With 84 direction bins, the 
ionisation front penetrates approximately 8 cells into the region 
shadowed by the clump, while with 42 direction bins this is al- 
most twice as many cells. The slice through the x-axis shows 
that the ionisation front outside the shadow is not affected by the 
dense clump. 

box size is 0.5/z _1 Mpc, the gas temperature is initially set to 100 
K. The sources are located in the 16 most massive haloes in the 
box, emitting f y = 250 ionising photons per atom over t s = 
3 Myr resulting in a photon flux of 



Ny=fy 



Qom H t s ' 



(37) 



where M is the total halo mass, Do = 0.27, = 0.043 and 
h = 0.7. The total simulation time is 0.4 Myr. SimpleX does 
not solve for the temperature state of the gas, so instead we as- 
sume a temperature of 1 • 10 4 K for the ionised gas. This might 
cause differences in the number of recombinations compared to 
the other codes in the comparison project that do solve for the 
temperature, since the recombination rate is a function of the 
temperature. 

5.3.1. Sampling function 

In order to perform this test, we first have to translate from the 
grid-based representation of the density field to the SimpleX grid. 
As discussed in Sect. I2.1.31 it is essential to have the highest 
dynamic range possible while keeping the density gradients to a 
minimum. Referring to Eq. (40) in KPCI09, we set the sampling 
parameter Q n to 5. We choose the parameter a in Eq. (TT31) close 
to its maximal value of 0.3 in order to have the highest resolution 
possible for this Q n in the low density regions. This results in 
no = 3.69 • 10" 5 cm~ 3 . The point density in the lowest density 
regions is equivalent to a resolution of approximately 77 3 for 
128 3 grid points. 

A slice through the z = Zbox/2 coordinate of the grid with the 
above parameters is shown on the right hand side in Fig. [TTJ If 
we compare the hybrid sampling scheme to the cubic sampling 




Fig. 18. Comparison between the result with the old sampling 
routine as was performed for the Radiative Transfer Comparison 
Project (left) and the hybrid sampling scheme (right). The num- 
ber of grid points is in both cases 128 3 , the mode of transport 
is ballistic only. For comparison reasons both grids have been 
interpolated to a regular grid of 128 3 cells. Shown is a slice 
through the z = Zbox/2 coordinate of the computational domain 
at t = 0.2 Myr, as the influence of the incorrect sampling is most 
pronounce there. In the left hand plot, the artefacts of under- 
sampling the low density regions are visible as neutral clumps, 
where no radiation has gone yet. The reason for this is the small 
number of grid points in these regions, resulting in a too small 
number of photons travelling there. The plot on the right hand 
side shows that this problem is solved by correctly sampling the 
low density regions. 



scheme, we can clearly see that the new scheme produces a grid 
that has the desired higher resolution in the dense filaments, but 
still has enough grid points in the low density regions to ensure 
that photons can travel into these regions. It is this grid that we 
will use for performing the radiative transfer simulations. 



5.3.2. The result of undersampling 

To stress the importance of sampling the medium correctly, we 
show in Fig. [18] the results using the old sampling method, 
used by SimpleX for the Radiative Transfer Comparison Project, 
and the hybrid sampling. For comparison purposes the mode of 
transport in both cases is ballistic and the number of grid points 
is 128 3 . The result of the incorrect sampling is clearly visible as 
dense neutral clumps in the ionised regions. This is not due to 
the fact that photons have preferential directions into the dense 
filaments, otherwise we would also see this effect in the hybrid 
sampling case. Rather, it is caused by the fact that there are too 
few grid points in the low density regions, resulting in very large 
cells. Since photons travel along the Delaunay edges, radiation 
simply does not travel into the low density regions, resulting in 
the observed large neutral clumps in the voids. 

An example of this effect is shown in Fig.[T9l Consider pho- 
tons travelling from left to right along the edges of the trian- 
gulation. The reason why the large cell is hard to ionise is not 
that there are not enough Delaunay edges pointing towards the 
cell. Clearly, the number of edges pointing towards the big cell 
is above average. However, as photons are sent to the d most 
straightforward neighbours and have no memory of their origi- 
nal direction, approximately half of the photons will be travel- 
ling around the big cell instead of ionising it. Thus, large cells 
are harder to ionise, resulting in the neutral clumps in the low 
density regions visible on the left-hand side of Fig. [TU This 
problem would be partly cured by the use of weights of the d 
most straightforward neighbours, giving the edges pointing into 
the big cell a higher number of photons. Another option would 
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Fig. 17. Left: Slice through the cosmological density field of llliev etaD ([2006) at coordinate z = Zbox/2. Centre: Cubic sampling 
using 128 3 grid points, as required by Eq. ([TIT) , showing that for realistic d ensity fluctuations the low and intermediate density 
regions are severely undersampled. Right: Sampling scheme described in Sect. 12.1.31 with parameters a = 0.3 and po = 3.69 • 10" 5 . 
The point density in the lowest density regions corresponds approximately to a resolution of 77 3 . 




Fig. 19. Example of the effect of large grid cells on radiation 
transport. Consider radiation travelling from left to right along 
the edges of the triangulation. Even though there are many 
Delaunay edges pointing towards the big cell, this cell will not 
be ionised easily. The reason for this is that radiation is sent to 
the d most straightforward neighbours and the photons have no 
memory of their original direction. In this case approximately 
half of the photons travelling from right to left that should be 
used to ionise the big cell, will instead be travelling around this 
cell. Weighting the most straightforward direction would allevi- 
ate this problem, as would restricting the opening angle in which 
the d most straightforward neighbours are allowed to be. 



be to restrict the opening angle in which the d most straightfor- 
ward neighbours are allowed to be, thus discarding most of the 
edges that point around the big cell in case of photons travel- 
ling from the right. Applying the direction conserving transport 
scheme will not solve this problem entirely, since the photons are 
still split up and travelling along the edges of the triangulation 
as with ballistic transport, so the problem is essentially the same. 
However, a slightly larger fraction of photons will be travelling 
into the big cell with direction conserving transport compared to 
ballistic transport, as photons that are send in a direction around 
the large cell remember their original direction and thus have 
a higher probability to travel into the direction of the big cell 
again. 



5.3.3. Grid dynamics 

We have performed a set of simulations to investigate the effect 
of regularly updating the grid according the changes in the op- 
tical depth. In Fig. [20| six different simulations are shown. We 
performed one simulation where no vertices were removed, one 
simulation where highly ionised vertices were removed until a 
minimum local resolution of 128 3 was reached and one sim- 
ulation where the minimum resolution of the vertex removal 
was 77 3 . Thus, the second simulation removes redundant opti- 
cally thin vertices until the resolution of the input hydrodynam- 
ics grid has been reached, removing the additional resolution of 
the adaptive grid when it's no longer necessary, while the third 
simulation removes redundant optically thin vertices until the 
minimum resolution of the SimpleX grid itself has been reached. 
In all three cases we also varied the number of direction bins in 
which the photons at optically thin vertices were stored between 
42 and 84 bins. 

When we compare the lower panels of Fig. |20l which show 
the position of the ionisation front at the end of the simulation 
time, a clear difference is visible between the simulations. The 
centre plot shows that the removal of redundant grid points un- 
til the resolution of the original hydrodynamics grid has been 
reached results in sharper shadows at the place of high density 
filaments, because the small amount of scattering due to the ran- 
dom rotations of the direction bins has been kept to a minimum. 
Moreover, in this simulation the difference between using 42 and 
84 direction bins is only a few percent, showing that the numer- 
ical scatter is indeed negligible. The effect of vertex removal is 
not apparent at earlier times in the simulation (figure [20j top 
panel). In this case, the ionised regions occupy a small volume 
in the computational box so the photons do not have to travel a 
large distance. Only at a later stage, when photons have to travel 
almost an entire box length, the influence of the numerical scat- 
ter is visible. The effect of the scatter remains small even at these 
late stages of the simulation, only visible near high density fila- 
ments. 

In the simulations shown on the right-hand side of Fig. [20| 
is visible that the minimum resolution should not be lower than 
the resolution of the original input grid. In that case we are not 
only removing the extra grid points that were put in the high 
density filaments due to the adaptive nature of the SimpleX grid, 
we are also removing grid points that are necessary to represent 
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Fig. 20. Test 3, slice through coordinate z = Zbox/2 of the cosmo- 
logical density field. Contours show the position of the ionisation 
front for simulations using 42 direction bins at optically thin ver- 
tices (white) and 84 direction bins (black). Shown are the results 
of simulations without vertex removal (top), vertex removal with 
a minimum local resolution of 128 3 (centre) and vertex removal 
with a minimum local resolution of 77 3 (bottom). The left panels 
show the ionisation front after 0.05 Myr, the right panels after 0.4 
Myr. The removal of vertices in highly ionised regions reduces 
the numerical scatter due to the random rotations of the direction 
bins, therefore the difference between the two simulations in the 
plot in the centre stays within a few percent. The centre plot also 
shows sharper shadows than the plot on the left-hand side. The 
right-hand side plot shows that when the minimum resolution is 
very low, the medium is not well represented anymore, causing 
the shadows to disappear. 









1 ' a 



Fig. 21. Slice through the z = Zbox/2 coordinate of the cosmolog- 
ical density field showing the H i fraction at time t = 0.05 Myr. 
Beginning in the top left hand corner and proceeding clock- wise 
are the results from C 2 RAY, crash, ftte and SimpleX. 




Fig. 22. Slice through the z = Zbox/2 coordinate of the cosmo- 
logical density field showing the H i fraction at time t = 0.4 Myr, 
the final simulation time. Beginning in the top left hand corner 
and proceeding clock- wise are the results from C 2 RAY, crash, 
ftte and SimpleX. Note that with the direction conserving trans- 
port, SimpleX is capable of producing sharp shadows behind 
dense structures. 



the physical structures in the original hydrodynamics grid. By 
removing too many grid points from the high density filaments, 
these structures are smeared out over a few large cells in the 
SimpleX grid. Hence, recombinations occurring in the filaments 
are spread out over too large a volume, resulting in a lack of 
shadowing in the dense filaments. 



5.3.4. Comparison to other codes 

The SimpleX run with vertex removal and a minimum local res- 
olution of 128 3 is compared to the results of other codes in the 
Radiative Transfer Comparison Project in Fig. |2T]and Fig. l22l 
Shown is the Hi fraction of a slice through the z = Zbox/2 coor- 
dinate of the computational volume at times t = 0.05 Myr and 
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t - 0.4 Myr. The res ults of the SimpleX sim ulation are compared 
to tho se of c 2 r ay dMellema et ai1l2006h. c rash (Ma selli et all 
120031) and ftte dRazoumov & Cardallll2QQ5|). Note that t he new 
version of the crash code, described in iMaselli et al.l d2009), 
shows results that are in better agreement with the C 2 RAY and 
ftte results. 

Both the position of the ionisation front as the inner structure 
of the ionised regions of the SimpleX simulation show very good 
agreement with the other codes, despite the fact that SimpleX 
does not solve for the temperature state of the gas. In the SimpleX 
simulation the temperature of the ionised gas is constant at 10 4 
K. This is slightly lower than the temperature that the other 
codes, which do solve for the temperature, find. Thus, SimpleX 
overestimates the number of recombinations that are occurring 
in the ionised gas. This results in an ionisation front that is 
slightly behind that of the other codes at the end of the simula- 
tion. Also visible in these figures is the effect of Poisson noise on 
the stochastic grid, causing a less smooth ionisation front. This 
effect decreases with higher resolution. We are currently investi- 
gating how to minimise the Poisson noise during the creation of 
the grid. 

With the direction conserving transport in the ionised re- 
gions, SimpleX is capable of reproducing the shadows behind 
neutral clumps and filaments that the other codes show as well. 
These shadows were absent in the simulation with ballistic trans- 
port only, the right hand side of Fig.|T8j On the whole, the mor- 
phology of the ionised region that SimpleX finds shows excellent 
agreement with the results of the other codes. 

6. Summary 

In this paper we have presented esse ntial updates to the orig - 
inal SimpleX algorithm described in iRitzerveld & Ickel (120061) . 
The most important improvement are the parallellisation for dis- 
tributed memory machines, the new sampling function and the 
transport mode for optically thin regions. 

All the radiation transport algorithms used in SimpleX are lo- 
cal, photons travel from one grid point to the other. This makes 
the method relatively straightforward to parallellise. We have 
shown that the parallel execution of SimpleX shows excellent 
scaling properties on distributed memory machines, the work 
load per processor is decreased by almost 50% when the num- 
ber of processors on which the algorithm is executed is doubled. 
An exception to this fortunate scaling is the triangulation algo- 
rithm, due to the extra points in the boundary around processors 
that need to be triangulated. For the present application this is 
not a problem, since the computation time of the triangulation is 
an order of magnitude smaller than the computation time of the 
radiative transfer itself. 

The new sampling functions allows us to effectively choose 
the gradient in point density as to get an idea of the numerical 
diffusion that might occur. Together with a choice for the mini- 
mum resolution in the low density regions, this results in a grid 
optimal for doing radiative transfer calculations with SimpleX. 
We have shown that by constructing the grid in the correct way, 
the artefacts that showed up with the original SimpleX algorithm 
when doing simulations of a realistic density field with large 
density fluctuations are absent. 

The direction conserving transport mode ensures that pho- 
tons keep their directions when they are travelling through re- 
gions with a low optical depth. This removes the numerical 
diffusion that occurs when photons are transported ballistically 
everywhere. Furthermore, with direction conserving transport 
SimpleX results show sharp shadows behind dense regions. Even 



though direction conserving transport is computationally more 
expensive than ballistic transport, by applying it only to the re- 
gions where it is necessary, the extra computational overhead is 
limited. 

We would like to stress that these improvements in the ac- 
curacy of the method have been achieved without sacrificing its 
main benefits. The algorithm remains computationally very ef- 
ficient, with a computation time independent of the number of 
sources. Due to the adaptive nature of the grid, the dynamic 
range of the resolution is very high. Combined with the paral- 
lellisation for distributed memory machines, simulations can be 
done with a very high numbers of grid points, making SimpleX 
an ideal tool for doing large scale reionisation simulations. 

We have shown that SimpleX produces very accurate results 
in various standard test problems for cosmological reionisation. 
In the classical case of a growing Hn region about a single 
source both the position of the ionisation front and the inner 
ionised region agree with the analytic expectations within 1%. 
In the test problem closest to our future application, that of a 
cosmological density field with multiple sources, the SimpleX 
results show an excellent agreement with the results of other 
codes. 
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